!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, po, x, y, nx, ny)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, po

   integer(4) :: i, j, i_total, j_total
   real(8) :: pr, dx1, dx2, dy1, dy2

!----------------------------------------------------------------------------------------------------------------------------------|
!  Shock dam
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx

!      if (abs (x(i) + 1.0d0) .lt. 0.025d0) then
      if ((t .gt. 0.00d0 .and. t .lt. 0.01d0) .or.                                                                                 &
          (t .gt. 0.20d0 .and. t .lt. 0.21d0) .or.                                                                                 &
          (t .gt. 0.40d0 .and. t .lt. 0.41d0)) then
         ro(i, j) = 1.5d0
         mx(i, j) = 1.0d0
         my(i, j) = 0.0d0
         mz(i, j) = 0.0d0
         bx(i, j) = 0.0d0
         by(i, j) = 0.0d0
         bz(i, j) = 0.0d0
         pr = 2.0d0
      else
         ro(i, j) = 1.0d0
         mx(i, j) = 0.0d0
         my(i, j) = 0.0d0
         mz(i, j) = 0.0d0
         bx(i, j) = 0.0d0
         by(i, j) = 0.0d0
         bz(i, j) = 0.0d0
         pr = 1.0d0
         if (x(i) .ge. 0.0d0 .and. x(i) .le. 0.4d0) then
            by(i, j) = (cos ((x(i) - 0.2d0) * pi / 0.2d0) + 1.0d0) / 2.0d0
            pr = 1.0d0 - by(i, j) * by(i, j) / 2.0d0
         endif
      endif

      en(i, j) = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                            &
         (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0 + pr / gmm1

      gx(i, j) = 0.0d0
      gy(i, j) = 0.0d0

      po(i, j) = 0.0d0

   enddo
   enddo

   if (t .eq. 0.0d0) then

      open (100, file = trim (directory) // '/density.dat', form = 'unformatted')
      read (100) nx_total_finest, ny_total_finest

      if (.not. allocated (x_total)) then
         allocate (x_total(nx_total_finest))
      endif

      if (.not. allocated (y_total)) then
         allocate (y_total(ny_total_finest))
      endif

      if (.not. allocated (ro_total)) then
         allocate (ro_total(nx_total_finest, ny_total_finest))
         read (100) x_total, y_total
         read (100) ro_total
      endif

      close (100)

      i = 1
      j = 1
      loop: do j_total = 1, ny_total_finest - 1
      do i_total = 1, nx_total_finest - 1
!         if (abs (x(i) - x_total(i_total)) .le. 1d-12 .and. abs (y(j) - y_total(j_total)) .le. 1d-12) then
         if (x(i) .gt. x_total(i_total) .and. x(i) .le. x_total(i_total + 1) .and.                                                 &
             y(j) .gt. y_total(j_total) .and. y(j) .le. y_total(j_total + 1)) then

            if (sum (ro_total(i_total : i_total + 1, j_total : j_total + 1)) .ge. 0.5d0 ) then
               dx1 = x(i) - x_total(i_total)
               dx2 = x_total(i_total + 1) - x(i)
               dy1 = y(j) - y_total(j_total)
               dy2 = y_total(j_total + 1) - y(j)

               ro(i, j) = (ro_total(i_total, j_total) * dx2 * dy2 + ro_total(i_total + 1, j_total) * dx1 * dy2 +                   &
                           ro_total(i_total, j_total + 1) * dx2 * dy1 + ro_total(i_total + 1, j_total + 1) * dx1 * dy1) /          &
                           ((dx1 + dx2) * (dy1 + dy2)) * 4.0d0 / 255.0d0 + 1.0d0

               pr = 1.0d0
               en(i, j) = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                      &
                  (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0 + pr / gmm1
            endif

            i = i + 1
            if (i .gt. nx) then
               i = 1
               j = j + 1
               if (j .gt. ny) exit loop
            endif
         endif
      enddo
      enddo loop

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
